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Abstract 

The second Bianchi identity can be recast as an evolution equation 
for the Riemann curvatures. Here we will report on such a system for 
a vacuum static spherically symmetric spacetime. This is the first of 
two papers. In the following paper we will extend the ideas developed 
here to general vacuum spacetimes. In this paper we will demonstrate 
our ideas on a Schwarzschild spacetime and give detailed numerical 
results. For suitable choices of lapse function we find that the system 
gives excellent results with long term stability. 

1 Introduction 

Despite a slow start, hyperbolic formulations of the Einstein equations have 
in recent times become the system of choice for numerical relativity. 

The confidence afforded to hyperbolic systems is borne out not just by the 
recent success in numerical relativity [1, 2, 3] but also from their strict math- 
ematical underpinnings (which guarantees that future evolutions exist and 
that they depend smoothly on the initial data, these are key aspects of the 
theorems that demonstrate the stability of the system, see [1] for details). 

One of the earlier hyperbolic formulations was given by Friedrich [■ >] in which 
he used the second Bianchi identities to evolve the Weyl curvatures in-situ 
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with the metric. This idea has been extended by many other authors [6, 7, 
8, 9] and the resuhing equations are commonly referred to as an Einstein- 
Bianchi system. 

Yet despite their mathematical elegance and the virtues that this would be- 
stow upon a numerical code there seems to be have been very few numerical 
applications employing an Einstein-Bianchi system (though see [10, 11]). 

In this paper we will report on a simple Einstein-Bianchi system adapted 
to a discrete lattice for static spherically symmetric spacetimes. We were 
lead to this formulation not by way of Friedrich's paper but rather as a 
direct extension of our own ideas developed in an earlier series of papers 
[12, 13, 14, 15]. In those papers we used the spatial form of the second Bianchi 
identities (i.e., the second Bianchi identity for the 3-metric) to compute the 
3-Riemann curvatures across a Cauchy surface. This device proved to be the 
key element in obtaining accurate and stable evolutions of the initial data. 

Our longer term intention is to employ an Einstein-Bianchi system to evolve 
a 3-dimensional lattice. This will require not only evolution equations for the 
legs of the lattice, such as those given in [13], but also evolution equations 
for the curvatures. This will be the subject of the second paper in this series. 

For the simple case presented here we find that the system works very well. 
The evolutions are stable, though this depends on the choice of the lapse 
function, see sections (4.1) and (6.1). We also find that the constraints are 
well behaved (they appear to grow linearly with time and converge to zero 
as the lattice is refined, see section (6.2)). 

As this paper borrows heavily from two of our previous papers, which we 
refer to as Paper 1 and Paper 2, we will skip over many of the derivations 
and arguments assuming instead that the reader is familiar with the material 
in Paper 1 and Paper 2. 

2 Spherically symmetric spacetimes 

In this paper we will be constructing lattice approximations to the Schwarzschild 
spacetime in various slicings. In each case the continuum metric can be writ- 
ten in the form 

ds"" = -N{r, tfdf + A{r, tfdr^ + S(r, tfd^^ (2.1) 
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for some set of functions N{r,t), A{r,t), B{r,t) and where dfl'^ = dd"^ + 
sin^ Odc/)'^ is the metric of the unit 2-sphere. We have introduced this coordi- 
nate form of the continuum metric simply as a precursor to the introduction 
of the lattice. As we shall soon see, we will use the coordinate lines and their 
local tangent vectors as a scaffold on which to build the lattice, after which 
we will have no further need for the coordinates (indeed we could dispense 
with the coordinates altogether at the possible expense of the clarity of the 
exposition) . 

Consider a local orthonormal tetrad built from the future pointing unit nor- 
mal to a typical Cauchy surface and three unit vectors m^, and 
where is parallel to the radial axis (see figure (2)). These basis vectors 
are also tangent vectors to the coordinate axes of a local Riemann normal 
frame. We will use this tetrad to record the frame components of the extrin- 
sic and Riemann curvatures on the lattice. Our notation, which we borrow 
from Paper 1, will be to use script characters to denote frame quantities, 
thus K-xx '■= Ki^.yTn'^my^ while TZtztz '■= R^lav|3't^^^"i^'^nz■ Also, to avoid an 
overflow of symbols, we will allow L^x and Lzz to represent both the length 
of the corresponding leg as well the leg itself. 

In this class of spacetimes, and on this tetrad, we know that the extrinsic 
curvature is diagonal and that a basis for the non-trivial Riemann curvatures 
is given by 

T^xyxyi T^xzxzi T^txtxi T^tztzi T^txxz 

Now using Rfj,i, = R'^^au and R = g^'^R^^ we find 

T^tz = —^IZtxxz (2-2) 

T^tt = T^tztz + 27^txtx (2-3) 

T^zz T^tztz ~l~ '^T^xzxz (2-4) 

T^xx T^txtx ~l~ T^xyxy ~l~ T^xzxz (2-5) 

= ~^T^txtx ~ '^T^tztz + 2 (TZxyxy + "^T^xzxz) (2-6) 

while the non-trivial vacuum Einstein equations yield 

TZtz = TZu = 7^,, = 7^,. = (2.7) 

Combining the above shows that we can express all of the non-trivial Rie- 
mann curvatures solely in terms of TZxyxy and TZxzxz, namely 

T^txxz = (2.8) 

T^tztz '^'T^xzxz (2'9) 

T^txtx T^xyxy ~\~ T^xzxz (2'10) 
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In obtaining these relations we used g'^'^ = —t^f^ + m'^m'^^ + mt^rriy + m'^m'^. 

Note that TZxyxy and TZxzxz are not independent for the simple equation _R = 
leads to 

T^xyxy ~\~ '^T^xzxz (^■-^-^) 

We will use this equation as a check on our numerical integrations (see section 
(6.2) for more details). 

3 The lattice 

The symmetries in the Schwarzschild spacetime allows us to use a very simple 
ladder-like structure for the lattice, as indicated in figure (1). One way to 
imagine the construction of the lattice is to consider the coordinate mesh 
generated by setting t = constant and 6 = 7r/2 in the coordinate form of 
the metric in (2.1). Then the rungs of the ladder are generated by small 
increments in leading to L^x ~ BA(j) while the side rails would coincide with 
two radial curves (i.e., = constant) with Lzz ~ AAr. Clearly, specifying all 
of the Lxx and L^z is equivalent to specifying the metric components A{r, t) 
and B{r,t). Note that throughout this paper we treat the Lxx and Lzz as 
continuous functions of time. 

We will label the nodes from to n and on the few occasions where we need 
to discuss more than one leg at a time we will write {Lxx)i to denote an Lxx 
leg at node i. In the same way {Lzz)i will denote the Lzz that joins the nodes 
i and z + 1. Similar notation will be used for other data on the lattice. 

The initial data (as described in section (5)) are constructed in a way that 
guarantees reflection symmetry at the throat (which is always tied to node 
0). 

In our computer code we extend our lattice a small way over the throat, by 
including the nodes -3 to -1, so that we can readily impose the reflection 
symmetries (by simply copying data across the throat, at no point do we 
independently evolve any of the data to the left of the throat). 

4 The evolution equations 

Our present task is to develop evolution equations for the leg-lengths, the 
extrinsic curvatures and, the principle innovation in this paper, evolution 
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equations for the Riemann curvatures. 

A simple derivation of the evolution equations for our lattice can be obtained 
from a general pair of equations developed in Paper 2. There it was shown 
that the first and second variations of arc lengths can be written in a form 
remarkably similar to the ADM equations, namely 

= -2ArK,.Axi;.Axr, + O (L^) (4.1) 



^ - 2Ar|,^Ax" Axf, (4.2) 



d f I dLl 

+ 2N {K^^K^p - R^a.^t^n Axf^. AxJ + O {L') 



Note that in the following we will ignore the leading error terms 0{L^). 
Applying these equations to the two legs L^x and Lzz of our spherically 
symmetric lattice leads immediately to 



dL^ 



dt 

dt 

dlQxQi 
dt 

dJCz, 
lit 



-NJCxxLxx (4.3) 

-NJCzzLz, (4.4) 

-Nxx + N {Utxtx + /CL) (4.5) 

-AT,,, + N [Uutz + /Cy (4.6) 



The last part of the picture is to provide evolution equations for the Riemann 
curvatures, TZxyxy and TZxzxz- The basic idea is to rearrange the terms in the 
Bianchi identities to isolate the time derivatives while estimating the spatial 
derivatives from data imported from neighbouring cells. The calculations are 
straight-forward but a bit tedious to present here so we defer the full details 
to the Appendix. This leads to the following evolution equations 

~ '2-^f^xx i^T^xyxy ~l~ T^xzxz) (4' 7) 



dt 

diT^xz. 
dt 



^xxT^xzxz ~\~ -^^zz {T^xyxy ~l~ '^T^xzxz^ (4-^ 



The Riemann curvatures TZxyxy and TZxzxz would normally not be evolved but 
rather derived from the lattice data such as the leg lengths Lxx and L^z- In 
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Paper 1 we used (discrete versions of) the geodesic deviation equation and 
the spatial Bianchi identity* 

= ^^^-/^^^^^^ - ^7^_,^ (4.10) 
dz dz 

to compute the 3- dimensional Riemann curvatures ^TZxyxy and ^TZxzxz on the 
lattice. In raising the TZxyxy and TZxzxz to dynamical variables on the lattice 
we are forced to view equations (4.9,4.10) as constraints on the lattice data. 
In section (6.2) we shall present discretised versions of these constraints which 
we will later use to check the quality of our numerical results. 

The one remaining constraint is the standard momentum constraint (see 
Paper 1 for details) 

diyLxx^xx) 1^ dLxx /. , -I \ 

^=—z ^'-''^ 



4.1 Artificial viscosity 

Our numerical experiments (which we will present shortly) showed that the 
future evolutions can be subject to high-frequency instabilities. This was 
seen to occur only in the cases where the lapse function was controlled by 
its own evolution equation (e.g., as in Harmonic slicing). For such cases 
we found that stability could be recovered with the addition of an artificial 
viscosity term to the evolution equations. 

Let W be any one of the dynamical variables, L^x-, L~,zi JC^x, ^zz- Then the 
artificial viscosity is introduced by the addition of a simple dissipation term 
to the evolution equation for W . After some experimentation we settled on 
the following form 

where dWi/dt is the right hand side of the original evolution equation (4.3- 
4.6) and /i is a (small) constant. Other choices were tried but this form 
seemed to produce stable evolutions for the longest periods of time. Note 
that we do not add the dissipation terms to the evolution equations for the 



*Here z is the proper distance measured up the middle of the lattice i.e., along a 
trajectory that passes through the mid-points of each La;^; 
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curvatures (doing so seemed to make no difference to the evolutions and had 
no effect in controlling the instabilities). 

How should fi be chosen? We need to choose it large enough to ensure 
that the evolution is stable over a given time interval while also keeping it 
sufficiently small so as to not effect the large scale features of the numerical 
solution. By trial and error we found that setting /i = 0.08 worked well for 
evolutions to t = 100.0 using n = 2048 nodes. We also found that as the 
number of nodes was increased we had to make a proportionate increase in 
fi to maintain the same quality of the evolution over the same time interval. 
That is fi = O (n). 

The dissipation term is easily seen to be a finite difference approximation to 
^Lzz<PW / dz^ and thus it may appear to be like a Kreiss-Oliger term that van- 
ishes in the continuum limit. However, since we are forced to set /i = (9 (n) 
and as Lzz = O [l/n) we see that the term ^Lzz is approximately constant, 
say /i', and thus the dissipation term is actually of the form n'dPW/dz^. This 
is a standard dissipation term commonly used in hydrodynamic simulations 
and it does not vanish in the continuum limit. 



5 Initial data 

The initial data on the lattice are the L^x, Lzz, ^xx, ^zz, T^xyxy and TZxzxz 
at each node of the lattice. Their time symmetric initial values were set by a 
combination of the Hamiltonian constraint, the geodesic deviation equation 
and the Bianchi identities. A full account of the choices made in coming 
to the equations described below can be found Paper 1. Here we will just 
quote the relevant equations simply to provide explicit details of how we 
constructed our initial data. 

To ensure that the initial data is time symmetric we set IC^x = and ICzz = 0. 

The Lzz were set according to the method of Bernstein, Hobill and Smarr [16] 
using n = 2048 on a grid of length 800.0m. The ADM mass, m, was set to 
be 1.0 and the L^x, Tlxyxy and TZxzxz, for i = 1, 2, 3 ... n, were set according 
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to 



{Lxx)i — {Lxx)i-1 + JJ \ {{Lxx)i-1 — {Lxx)i-2) 
\J^zz)i-2 



— -^{Lzz)i-1 {{Lzz)i-1 + (-^22)1-2) {LxxT^xzxz)i^i (^-l) 



('^a;2a;2)i — {T^xzxz)i-1 ( r^T2 \ (T2 \ ^ (^-2) 

\<3 \^xx)i ~ K^xx)i-l/ 

(T^xyxy^i '^{^'^xzxz) i (^•'^) 

At the reflection symmetric throat (i.e., at node 0) we set {Lxx)o — 0.01 and 

{'^xyxy)o — ~'^{T^xzxz)o = 0.25. 



6 Results 

In all of our results we used a 4-th order Runge-Kutta integrator with the 
time step set equal to 1/2 the smallest L^z on the lattice (which happens to 
be {Lzz)o)- 



6.1 Slicing conditions 

We ran our code for eight distinct slicing conditions, some were set by simple 
algebraic expressions while others involved differential operators. 

We made four choices for the algebraic shcings, 

= exp(-2/C^^) (6.1) 

N = (6.2) 

N = exp{-nxyxy) (6.3) 

N = z ^ (6.4) 

^ ' '^xyxy 
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and three choice for the differential shcings, 

dN 

1 + log shcing — = -2NK (6.5) 

at 

Harmonic slicing — — = — K (6.6) 

(Jjb 

Maximal shcing V^iV = ^RN (6.7) 

The eighth slicing condition was the simple case of geodesic slicing = 1. 

The algebraic slicings were introduced after our early explorations with the 
differential lapses, all of which developed high-frequency instabilities after 
a short time (well before t = 100). The algebraic slicings did not require 
any artificial viscosity and performed remarkably well, showing no signs of 
instabilities to at least t = 1000 (excluding the lapse (6.2) which hits the 
singularity at t ~ 32). We have not run our codes beyond t = 1000 so we 
can not comment its stability for t > 1000. 



6.2 Code tests and results 

We subjected our code to many of the tests used in Paper 1, such as the time 
at which geodesic slicing hits the singularity, the rate at which the lapse at 
the throat collapses in maximal slicing and the constancy of L^x on the hori- 
zon. The results for these various slicings are shown in figures (4-9). All of 
the results are as expected. For the geodesic slicing the code crashes at ap- 
proximately one time step short of the singularity. The familiar exponential 
collapse of the lapse for maximal slicing is evident in figure (6). In this case 
it is known that the lapse at the throat should behave as A^ ~ (3exp{at) for 
t oo with a = -(2/3)(^/2) ~ -0.54433, see [17]. We estimated the slope 
of In A^ vs t from our numerical data to be —0.54215 which agrees with the 
exact value to within 0.4 percent. 

We also have a new test obtained by a simple combination of the evolution 
equations. From equations (4.7,4.8) we find that 

xzxz ) (6.8) 

and as = TZxyxy + '^T^xzxz on the initial slice (by construction, see (5.3)) we 
conclude that = TZxyxy + '^T^xzxz for all time. This is not surprising, our 
evolution equations for the curvatures are based on the Bianchi identities 
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and these are guaranteed to preserve the constraints. If we now set 
T^xyxy + "^T^xzxz (4-8) and combine the resuh with (4.3) we find 



^^xx i^xyxy 

dt 



(6.9) 



This gives us a new test of our code, that the quantity Ll^TZ^yxy should 
be constant throughout the evolution. Importantly this applies to all slicing 
conditions. In figure (9) we have plotted the fractional variations in L^^TZxyxy 
for two choices of slicings. We see that the errors for the 1 + log slicing are 
much larger than those for the algebraic slicing which we attribute to the 
use of an artificial viscosity. This last claim is easily checked by varying the 
artificial viscosity parameter /z. We find that the errors in L^^TZ^y^y varies 
linearly with fi. Note that in obtaining equation (6.9) we have ignored the 
higher order error terms that would arise if we had carried through the O {L^) 
truncation error from (4.1). Thus even if we set /i = we can expect some 
variation of L^^TZ^yxy over time (though this variation should vanish more 
rapidly than O {L^)). 

We also have three constraint equations, namely the geodesic deviation equa- 
tion (4.9), the 3-dimensional Bianchi identity (4.10) and the momentum con- 
straint (4.11). The discrete form of these equations are 

Q = ^ """"^ - 'TZxzxz^^ (6.11) 
^ Dz Dz ^ ' 

Dz Dz ^ ' 

where TZ^zxz is the average of TZxzxz across L^z while D/Dz and D/Dz are 
discrete derivative operators defined as follows. For a typical smooth function 
f{z) sampled at the grid points Zi we define 



(6.13) 




J i ^zz + ^zz V V J V ^zz 

D'f\ 2 fr-r r-f 



(6.15) 
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where we have introduced the superscripts +, o and - to denote quantities 
at the grid points Zj+i, Zi and Zi_i respectively. Note that the sample points 
Zi are constructed from the lattice L^^ by the recurrence relation Zi^i = 
Zi + {Lzz)i with zq = 0. In this notation we have TZ^yxy '■= {T^xyxy^^xyxy) 1'^- 
Finally we note that the 3-curvatures can be computed from the 4-curvatures 
by way of the Gauss equation, 

T^xyxy = T^xyxy ~ ^xx (6.16) 
T^xzxz T^xzxz ^xx^zz (S'l'^) 

Ideally we would like to see P = Q = M = but in reality we expect Pi, 
Qi and Mi to be non-zero but small. This is indeed what we observe, see 
figure (9). We also computed a crude estimate of the rate of convergence (of 
Q, P and M to zero at a fixed time) by running our code twice, once with 
n = 2048 and once with n = 1024 and then forming suitable ratios of the 
constraints at the horizon. In this manner we estimated, in the absence of 
artificial viscosity, that P = O {n~^), Q = O [n^"^) and M = O {n~^) while 
the addition of artificial viscosity degraded the convergence to P = O {n~^), 
Q = (n-^) and M = C {n''^). 

We also tried setting D/Dz := D/Dz and TZ^yxy '■= TZ-xyxy discrete 
Bianchi constraint but this lead to a reduction in the rate of convergence. 
The form of the discrete Bianchi constraint as given above (6.11) is readily 
seen [I I] to be a second-order accurate estimate to the continuum Bianchi 
identity at the centre of the leg L^z- 

One might ask why we have not included the Hamiltonian constraint in 
our code tests. The simple answer is that it is trivially satisfied by our 
discrete equations. This follows from the discussion surrounding equation 
(6.8) where we showed that = TZxyxy + '^TZxzxz for all time. It follows that 
the Hamiltonian H := G^ytH^ will also vanish for all time. Note that this 
analysis was based on our discrete equations, not on the continuum equations. 
We did indeed check that our code maintained = 'R-xyxy+'^Tlxzxz throughout 
the evolution. 
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A Bianchi identities 



Here we will use the Bianchi identities to obtain evolution equations for the 
two curvatures TZxyxy and TZxzxz- We will follow the method given in our 
earlier paper [1 1] in which we used data imported from the neighbouring 
computational cells to estimate (by a finite difference approximation) the 
various derivatives required in the Bianchi identities. We will employ Rie- 
mann normal coordinates^, one for each computational cell, with the origin 
centred on the central vertex and the coordinate axes aligned with those de- 
scribed in section (2), see also figure (2). In these coordinates, the metric in 
a typical computational cell is given by 



where L is a typical length scale for the computational cell and g^^^, and 
R^auii are constant throughout the computational cell. A convenient choice 
for g^y is diag(— 1, 1, 1, 1) (such a choice can always be made by suitable gauge 
transformations within the class of Riemann normal frames). In this case the 
frame components TZxyxy and TZxzxz reduce to the coordinate components 
Rxyxy and Rxzxz rcspcctively. A further advantage of using Riemann normal 
coordinates is that at the origin, where the connection vanishes, covariant 
derivatives reduce to partial derivatives. 

The two Bianchi identities that we need are 



This pair of equations contains 4 spatial derivatives each of which we will 
estimate by a finite difference approximation. But in order to do so we must 
first have a sampling of the 4 curvatures at a cluster of points near and 
around the central vertex. Our simple ladder-like lattice, with its collection 
of computational cells along one radial axis, would allow us to compute only 
the z partial derivatives. For the x and y derivatives we will need to extend 
the lattice along the x and y axes. In short we need a truly 3 dimensional 
lattice. Fortunately this is rather easy to do for this spacetime. We can use 
the spherical symmetry of the Schwarzschild spacetime to clone copies of the 
ladder (by spherical rotations) so that a typical central vertex of the parent 

^^For more details on Riemann normal coordinate see [18] and the references cited 
therein. 




X 



x^x'< + O (L^) 





(A.2) 
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ladder-lattice becomes surrounded by 4 copies of itself. It has two further 
nearby vertices, fore and aft along the radial axis, that are themselves central 
vertices of neighbouring cells in the original ladder-like lattice. In figure (3) 
we display an xz slice of the cloned lattice. 

We now need the coordinates of all six of the neighbouring vertices. This 
would require a solution of 

Ll = {x^ - (x^^ - xr) - + O {L') (A.3) 

for the xf for given values for the Lij and R^avis- However, as we are only 
going to use these coordinates to construct transformation matrices which 
will in turn multiply the Riemann curvatures, it is sufficient to solve (A.3) 
using a fiat metric. Note that the above equations can only be used to 
compute (in fact estimate) the spatial coordinates of the vertices. For the 
time coordinates we can appeal to the smoothness of the underlying metric"'' 
to argue that for each vertex t = O {L"^)- The result is that the typical central 
vertex, with coordinates (0, 0, 0, 0), will have 6 neighbouring central vertices 
with coordinates as per Table (1). 



Vertex 


t 




y 


z 





( , 


, 


, 


) 


1 


( , 




, 


) 


2 


( , 


, 


Lyy , 


) 


3 


( , 


Lixx 1 
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) 


4 


( , 


, 


~Lyy, 


) 


5 


( , 


, 


, 


Lzz ) 


6 


( , 


, 


, 





Table 1: Riemann normal coordinates, to O {L'^), of the central vertex and its 
6 immediate neighbours. These coordinates were computed using a flat space 
approximation. 

This accounts for the structure of our lattice but what values should we 
assign to the curvatures at the newly created vertices? Let (A) pq denote the 

■^If (t, x*) are the coordinates for a local Riemann normal frame, then a smooth Cauchy 
surface through (0,0,0,0) is described locally by 2t = —KijX^x^ and as each x"^ — O (L) 
we also have t = O (L^) . 



13 



value of a quantity A at the vertex P in the local Riemann normal frame for 
vertex Q. Since our spacetime is spherically symmetric we can assert that 



{A)oo = {A)n = {A)22 = {Ahs = (A) 
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Then the idea that we will import data from neighbouring cells can be ex- 
pressed as 

iA)pQ = iU)pQ{App) 

where (U) pq is the transformation matrix, evaluated at P, from the Riemann 
normal frame of P to that of Q. This matrix will be composed of spatial 
rotations and boosts. 

To get the correct estimates for the first partial derivatives we need only 
compute U to terms linear in the leg-lengths. 

As an example, let us suppose wc wished to compute v^^^ for a spherically 
symmetric vector field v on the lattice. We start with (t')io = (^)io('^)ii and 

(C/)io = (S)io(i?)io 

where (-R)io represents a rotation in the x — y plane and (-B)io a boost in the 
t — X plane. Note that as we are working only to linear terms in the lattice 
scale the order in which we perform the rotation and boost does not matter. 
Thus we have 



(R) 



10 





COS a sin a 











— sin a cos a 



cosh /3 sinh ,5 

sinh/3 cosh/3 

10 

1 



The columns in the above matrices are labelled {t,x,y.z) from left to right 
and likewise for the rows. As we will latter be forming products of these 
matrices with the curvatures it is sufficient to compute these matrices as 
if we were working in flat spacetime. Thus to leading order in the lattice 
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spacing we find^ 



cosa = 1 + C (L^) , sin 



a 



dz 



cosh /3 = 1 + C (L^) , sinh /3 = -K^^L^^ + O (L^) 



and thus 



1 

- K T 





In a similar manner we find 



20 



1 KyyLyy 



^yy-^yy 




(f/)30 = (i?)30(i?)30 



XxL XX 













dLxx 


1 






dz 





1 





dLxx 
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dz 





K T 

^^xx-'-'xx 
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dLyy 




dz 





dLyy 


1 


dz 




KxxLxx 







- 


dLxx 


1 






dz 





1 





dLxx 









1 


dz 





+ O {L') 



(f/)40 = (-B)4o(-R)40 



1 









KyyLyy 

1 



^yy-^yy ^ 







1 

dL 



yy 



dz 






dL 



yy 



dz 
1 



+ O {L') 



§For the rotations we use standard Euclidian trigonometry, for the boost we use the 
definition n'^ — = —K^^{x\ ~ Xj) where is the future pointing unit normal to the 
Cauchy surface at the point a. 
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For the remaining two matrices, (t/)5o and {U)go, the job is quite simple, 
these matrices are built solely on boosts. This leads to 



(^7)50 = {BMR) 



50 



1 







1 

1 





'KzzLz 



1 



{U)m — {B)eo{R)eo 



1 




K,,L,, 



1 



1 









LT T 
^^zz-'-'zz 







+ O {L') 



Returning now to the construction of (f )io, we have 

K)io - {u^.hoivnii 



J^xx^xxV , —Kxx^xxV H ; V , U, ; V 

dz dz 
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and 



^^xx^xx'-' -I ^^xx^xx'^ J ^ ) - ^ 



dz dz 

We are now in a position to finally compute (t'*r)oO) to wit 



33 



00 



30 



2L, 



(t^*)ll - (t^033 ^ K)n + (t^-)33 , .^..g. 

Here we have written the truncation errors as O (L") with a > for it is not 
clear, at this level of analysis, what the exact nature of this term is (save 
that it vanishes as L — > 0). Since our spacetime is spherically symmetric we 
have 

{v)oo = (^^)ll = {v)22 = ('y)33 = {v)44 
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and thus 

Similar calculations can be used to compute all of the spatial derivatives of 
at the central vertex. 

We can now return to the principle objective of this section - to compute the 
various partial derivatives of the curvatures. We proceed exactly as above 
but with a minor change in that we will no longer carry the truncation errors 
within the calculations. Thus we have 

{Rixav^)i0 — {Ui/)io{Ua'')io{Ui,^)io{Ui3^)io{RT-pS\)ii 

for 1 = 1,2, 3,4, 5, 6 and {U^%o = g^.^g^'^iU^U with g^, = diag(-l, 1, 1, 1). 
And, as before, 

due to spherical symmetry. 

Using the above expressions for the {U)io and the following finite difference 
approximations 

/D \ _ i^tyxy) 10 - {Rtyxy) 30 

\-'^tyxy,x)00 — 



{Rtxxy,y)oO — 
{Rtzxz,x)oO — 
{Rtxxz,z)oO — 





9 r 

^■'-'xx 


{Rtxxy , 


)20 — {Rtxxy)40 




9 T 

yy 


{■^tzxz, 


ho — (Rtzxzjso 




2Lxx 


{■^txxz, 


)50 ~ {Rtxxz)m 



2L 



zz 



we find that 



■Rtyxy,x ^xx {^xyxy ~l~ Rtyty) ~l~ ^ 7 Rtyyz (-^•^) 



Rtxxy,y ^yy {Rxyxy ~l~ Rtxtx) j , Rtxxz (-''^•^) 



Lxx 


dz 


1 


dLyy 


^yy 


dz 


1 


dLxx 


Lxx 


dz 



Rtzxz,x ^xx {Rxzxz ~l~ Rtztz) ~l~ T 7 Rtxxz (-^•^) 

XX t///^ 

Rtxxz,z — ^zz {Rxzxz ~l~ Rtxtx) {■^•'^) 
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We have dropped the 00 subscript as we no longer need to distinguish between 
the neighbouring frames. By spherical symmetry we have 



Lxx -^yy 1 -^xx -^yy i -Rtxxz Rtyyz 1 Rxzxz Ryzyz 

while from the vacuum Einstein equations we have 






— Rtz — 


Rtxxz Rtyyz 







Rxx 


Rxyxy ~l~ Rxzxz 


Rtxtx 





= Ryy = 


Rxyxy ~\~ Ryzyz 


Rtyty 





Rzz 


Rxzxz ~\~ Ryzyz 


Rtztz 



Combining the last few equations leads to 

Rtxtx Rtyty Rxyxy ~l~ Rxzxz 
Rtxxz Rtyyz , Rtztz '^Rxzxz 

Substituting these into the above equations (A.4-A.7) and subsequently into 
the previous expressions for the Bianchi identities (A.1,A.2) leads to the 
following pair of equations 

dRxyxy _ r)T^ fr)T3 I p N 

- z,l\xx y^J^xyxy ~r J^xzxz) 



dt 
dR 

xz: 

dt 



'^KxxRxzxz ~l~ Kzz i^Rxyxy ~l~ '^Rx 



Our job is almost complete, but we still have two tasks ahead of us i) to 
introduce a lapse function and ii) to account for the limited time interval 
over which a single Riemann normal frame can be used. The first task is 
rather easy, we simply make the coordinate substitution t — )■ Nt leading to 

- INK XX ^Rxyxy ~l~ Rxzxz) (-^•^) 



dt 
dt 



3N KxxRxzxz ~l~ K zz i^Rxyxy ~l~ '^Rxzxz) (-^■^) 



and where we now have {gij,u)o = diag(— A^^, 1, 1, 1). The lapse A^ can be freely 
chosen at each vertex of the lattice (but subject to the obvious constraint 
that A^ > 0). The second task is a bit more involved. We know that each 
Riemann normal frame is limited in both space and time. Thus no single 
Riemann normal frame can be used to track the evolution for an extended 
period of time. We will have no choice but to jump periodically to a new 



18 



frame. This can be elegantly handled in the moving frame formalism. Thus 
our task reduces to finding a new set of evolution equations for the frame 
components TZxyxy and TZxzxz based on the equations given above for Rxyxy 
and Rxzxz- 

Let e^ai a = t,x,y, z be an orthonormal tetrad^, tied to the worldline of the 
central vertex and aligned to the coordinate axes. Thus we have e'^t as the 
future pointing tangent vector to the worldline while e'^z points along the 
z-axis. Then 

dt ~c/t {RfiavpG^xG yG x^^ y) 
(it c/t (-^/ja!^/9^'^a;^ 2^ x^^ 

Since our spacetime is spherically symmetric it is not hard to see that the 
tetrads of two consecutive cells (on the vertex worldline) are related by a 
boost in the t — z plane (arising from gradients in the lapse function). A 
simple calculation shows that 

dt ^ dt ' dt ,z z , ,z 
which when combined with the above leads to 



dT^xyxy 1 


( dR^OLUP 


dt ~ ' 


y dt 


dT^XZXZ i 


( dR^OLUP 


dt ~ ' 


[ dt 



6 y 

(A.IO) 

^ x^ 26 x^ z 2 ^ R^aul3(^^t(^ x^ 2 (-^■-^-^) 

In our frame we have chosen {g^u)o = diag(— A^^, 1, 1, 1), e'^a = ^ — 

X, y, z and e^t = l/A^, thus we see that the last term in the previous equation 
is proportional to TZtxxz- But for the Schwarzschild spacetime we know that 
T^txxz = and thus we have 

x^ y^ x^ y (■^■-^2) 



dl^xyxy i 


f dR^aul3 


dt ~ ' 


\ dt 


dT^XZXZ I 


( dR^auP 


dt ~ ' 


[ dt 



x^"^ z^^ x^^ z (A.13) 



which, when combined with (A. 8, A. 9), leads immediately to the evolution 
equations (4.7,4.8) quoted in section (4). 

^This tetrad is identical to that used in section (2), the change of notation introduced 
here is simply to avoid unwanted clutter in the following equations. 
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Figure 1: A simple lattice for a Scliwarzschild spacetime. This consists of two 
identical halves joined at the throat (denoted by the thick line). In our computer 
code we only store the right hand half (plus a few nodes from the left half to ensure 
reflection symmetry at the throat). 



y 




Figure 2: The local structure of the lattice. The (x, y, z) are a set of coordinates 
local to this set of legs. There is one such coordinate frame for each Lxx along the 
lattice. These coordinates are never used in the computer code but help to define 
the metric in the neighbourhood of L^x- 
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Figure 3: This shows an xz section of the lattice obtained by cloning the original 
2-dimensional lattice. The data in the upper and lower cells are identical to that 
in the middle cell, this follows from spherical symmetry. The small squares denote 
the central vertices of each computational cell. The angle a can be computed using 
standard Euclidean trigonometry as described in the text. 
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2 . OOe+00 



1.50e+00 - 



1. OOe+00 - 



5.00e-01 



. OOe+00 




1.0 1.5 2.0 2.5 

Proper time r 



3.0 3.5 




-4.00e-10 



Figure 4: The Schwarzschild areal coordinate r and the lattice L^x at the throat 
are related by r{t) = 2mLxx{t)/Lxx{0). In geodesic slicing r(t) at the throat is 
described by the parametric equations r{t) = m(l + cos ri{t)), t{ri) = m{r] + sinr/). 
These equations allow us to plot the exact evolution of r{t) (the smooth curve) 
against estimates from the lattice (solid points). The relative errors arc seen to be 
very small and are dominated by the truncation errors in the Runge-Kutta scheme. 



22 




Figure 5: The curvature profiles for four choices of lapse function. Each figure 
shows the lapse function for t = to t = 100 in steps of 10. The small diamond 
on each curve represents the location of the apparent horizon. There are only two 
curves visible in the algebraic slicing N = 20Lxx/ (1 + '^OLxx) due to the rapid rise 
in the curvatures as the slicing approaches the singularity. 
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Maximal lapse 



l.OOe-80 
l.OOe-90 




1 + log lapse 

dN/df = -2NIC 



20.0 40.0 60.0 

Radial proper distance 



100-0 



20.0 40.0 60.0 80.0 

Radial proper distance 




Figure 6: As per figure (5) but this time we display the lapse profiles. Note how 
quickly the lapse collapse at the throat in the 1 + log slicing. This would likely 
cause serious underflow problems for t > 300. Notice also the uniform spacing 
of the curves along the logarithmic axis for the maximal lapse. This show that 
the lapse collapses exponentially at the throat (a well known result for maximal 
slicing) . 
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Figure 7: The Riemann curvature TZxyxy and the lapse function for the alge- 
braic slicing with = 1/(1 + TZxyxy)- The left pair of figures shows the evolution 
over the range t = to 100 in steps of 10, while the right pair cover the range t = 
to 1000 in steps of 100. This lapse does not appear (on this time scale) to exhibit 
an exponential collapse at the throat. The profiles for the algebraic slicings appear 
to propagate into the asymptotically flat regions far more rapidly than any of the 
differential slicings. This may be explain why an artificial viscosity was not needed 
for the algebraic slicings - they carry away any small numerical errors before they 
have chance to grow. 
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Figure 8: In the absence of gravitational radiation the area of the horizon should 
remain constant. It follows that the Lxx should be constant on the horizon. Here 
we plot the fractional variation of L^x on the horizon. The irregular behaviour of 
the plots for later times is due, in part, to the difficulty in accurately locating the 
horizon. 
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Figure 9: According to equation (6.9) the quantity L^^TZxyxy should be conserved. 
In the top row we display the relative error, defined by 1 — C(t) /C(0) with C{t) = 
Lxx{t)Tlxyxy{t), for two choiccs of lapse. The errors for the 1 + log slicing are much 
larger than those for the algebraic slicing which we attribute to the action of the 
artificial viscosity terms. In the bottom row we display the momentum constraint 
for the same pair of slicings. This shows a slow growth in the momentum constraint 
over time (judging by the peaks in the plots, the growth appears to be linear in 
time). 
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